MODULE fefmm_LIP2_Fij

   USE fmm_global_paras
   IMPLICIT NONE
   PRIVATE
   ! Public procedures
   PUBLIC :: fefmm_build_Fij_LIP2

!   REAL(REALK),   PARAMETER :: T_VALUE_THR = 0.5D0
   INTEGER(INTK), PARAMETER :: LIPN = 2

CONTAINS

!-------------------------------------------------------------------------------

   SUBROUTINE fefmm_build_Fij_LIP2(FEdim,NFS,FEgrain,t,Fij)

      IMPLICIT NONE
      INTEGER(INTK), INTENT(IN)  :: FEdim, NFS
      REAL(REALK),   INTENT(IN)  :: FEgrain, t
      REAL(REALK),   INTENT(OUT) :: Fij(:,:)

      INTEGER(INTK) :: i,k, u,v

      Fij = zero

      DO i = 1, FEdim 
         v = i + FEdim*((NFS-1)/2)  ! V-point

         DO k = 1, NFS*( (FEdim -1)/LIPN )   ! finite elements
            u = 1 + LIPN*(k-1)    ! basic rho-point

            ! 1st LIP
            Fij(u,i)   = Fij(u,i)   + get_F1_LIP2(t,FEgrain,v,u)
            ! 2nd LIP
            Fij(u+1,i) = Fij(u+1,i) + get_F2_LIP2(t,FEgrain,v,u)
            ! 3rd LIP
            Fij(u+2,i) = Fij(u+2,i) + get_F3_LIP2(t,FEgrain,v,u)

         END DO

      END DO

   END SUBROUTINE fefmm_build_Fij_LIP2

!-------------------------------------------------------------------------------

   FUNCTION get_F1_LIP2(t,lambda,v,u)

      IMPLICIT NONE
      REAL(REALK),   INTENT(IN) :: t
      REAL(REALK),   INTENT(IN) :: lambda
      INTEGER(INTK), INTENT(IN) :: v  ! this is the V-point
      INTEGER(INTK), INTENT(IN) :: u  ! this is the rho-point
      REAL(REALK) :: get_F1_LIP2

      INTEGER(INTK) :: k
      REAL(REALK)   :: F1, F2

      k = v-u

      F1 = (exp(-t ** 2 * k ** 2 * lambda ** 2) * k - 0.3D1 * exp(-t ** 2    &
         * k ** 2 * lambda ** 2) + exp(-t ** 2 * lambda ** 2                 &
         * (k - 0.2D1) ** 2) - exp(-t ** 2 * lambda ** 2 * (k - 0.2D1) ** 2) &
         * k) / lambda / t ** 2 / 0.4D1

      F2 = (((-k ** 2 / 0.4D1 + 0.3D1 / 0.4D1 * k - 0.1D1 / 0.2D1)           &
         * erf(-dble(2 * t * lambda) + dble(t) * k * dble(lambda))           &
         + (k ** 2 / 0.4D1 - 0.3D1 / 0.4D1 * k + 0.1D1 / 0.2D1)              &
         * erf(dble(t) * k * dble(lambda))) / dble(t) + (0.1D1               &
         / dble(lambda ** 2) * erf(dble(t) * k * dble(lambda)) / 0.8D1       &
         - 0.1D1 / dble(lambda ** 2) * erf(-dble (2 * t * lambda) + dble(t)  &
         * k * dble(lambda)) / 0.8D1) / dble(t ** 3)) * sqrt(0.3141592654D1)

      get_F1_LIP2 = F1 + F2

   END FUNCTION get_F1_LIP2

!-------------------------------------------------------------------------------

   FUNCTION get_F2_LIP2(t,lambda,v,u)

      IMPLICIT NONE
      REAL(REALK),   INTENT(IN) :: t
      REAL(REALK),   INTENT(IN) :: lambda
      INTEGER(INTK), INTENT(IN) :: v  ! this is the V-point
      INTEGER(INTK), INTENT(IN) :: u  ! this is the rho-point
      REAL(REALK) :: get_F2_LIP2

      INTEGER(INTK) :: k
      REAL(REALK)   :: F1, F2

      k = v-u

      F1 = (-k * exp(-t ** 2 * k ** 2 * lambda ** 2) + 0.2D1               &
         * exp(-t ** 2 * k ** 2 * lambda ** 2) + k * exp(-lambda ** 2      &
         * t ** 2 * (-0.2D1 + k) ** 2)) / t ** 2 / lambda / 0.2D1

      F2 = (((-k + k ** 2 / 0.2D1) * erf(-dble(2 * lambda * t)             &
         + dble(t) * k * dble(lambda)) + (k - k ** 2 / 0.2D1)              &
         * erf(dble(t) * k * dble(lambda))) / dble(t)                      &
         + (-0.1D1 / dble(lambda ** 2) * erf(dble(t) * k                   &
         * dble(lambda)) / 0.4D1 + 0.1D1 / dble(lambda ** 2)               &
         * erf(-dble(2 * lambda * t) + dble(t) * k                         &
         * dble(lambda)) / 0.4D1) / dble(t ** 3)) * sqrt(0.3141592654D1)

      get_F2_LIP2 = F1 + F2

   END FUNCTION get_F2_LIP2

!-------------------------------------------------------------------------------

   FUNCTION get_F3_LIP2(t,lambda,v,u)

      IMPLICIT NONE
      REAL(REALK),   INTENT(IN) :: t
      REAL(REALK),   INTENT(IN) :: lambda
      INTEGER(INTK), INTENT(IN) :: v  ! this is the V-point
      INTEGER(INTK), INTENT(IN) :: u  ! this is the rho-point
      REAL(REALK) :: get_F3_LIP2

      INTEGER(INTK) :: k
      REAL(REALK)   :: F1, F2

      k = v-u

      F1 = (k * exp(-t ** 2 * k ** 2 * lambda ** 2) - exp(-t ** 2         &
         * k ** 2 * lambda ** 2) - exp(-lambda ** 2 * t ** 2              &
         * (-0.2D1 + k) ** 2) - k * exp(-lambda ** 2 * t ** 2             &
         * (-0.2D1 + k) ** 2)) / t ** 2 / lambda / 0.4D1

      F2 = (((k / 0.4D1 - k ** 2 / 0.4D1) * erf(-dble(2 * lambda * t)        &
         + dble(t) * k * dble(lambda)) + (-k / 0.4D1 + k ** 2 / 0.4D1)       &
         * erf(dble(t) * k * dble(lambda))) / dble(t)                        &
         + (0.1D1 / dble(lambda ** 2) * erf(dble(t) * k                      &
         * dble(lambda)) / 0.8D1 - 0.1D1 / dble(lambda ** 2)                 &
         * erf(-dble(2 * lambda * t) + dble(t) * k * dble(lambda)) / 0.8D1)  &
         / dble(t ** 3)) * sqrt(0.3141592654D1)

      get_F3_LIP2 = F1 + F2

   END FUNCTION get_F3_LIP2

!-------------------------------------------------------------------------------

END MODULE fefmm_LIP2_Fij
